//  only for private communication
//
//  Xianguo Lu 
//  lu@physi.uni-heidelberg.de
//  Xianguo.Lu@cern.ch
//

#ifndef XGLUTILS_H
#define XGLUTILS_H

#ifndef TMATRIXDSYM_H
#include "TMatrixDSym.h"
#endif

class TAxis;
class TCanvas;
class TChain;
class TF1;
class TGraph;
class TGraphAsymmErrors;
class TGraphErrors;
class TH1;
class TH1D;
class TH2D;
class TLatex;
class TLegend;
class TMinuit;
class TPad;
class TTree;

class XGLUtils
{
 public:
  enum{
    KPNG = 1,
    KPDF = 2,
    KEPS = 4,
    KROOT = 8,
    KC = 16,
    KJPG = 32,
    KPS = 64,
    KNORX = 0,
    KNORY = 1
  };

 private:
  XGLUtils(const XGLUtils &p);
  XGLUtils & operator=(const XGLUtils &p);

  //===================================================================
  //                                Functions
  //===================================================================
 public:
  typedef Double_t (*FFunc)(const Double_t *xx, const Double_t *par);

  static Double_t ProtectedExp(const Double_t xx);
  static Double_t ProtectedPow(const Double_t xx, const Double_t yy);

  static Double_t QuadSum(const Double_t x1, const Double_t x2);
  static Double_t QuadSum(const Double_t x1, const Double_t x2, const Double_t x3);

  static Double_t InvSqrt(const Double_t *xx, const Double_t *par);
  static Double_t PolN(const Double_t *xx, const Double_t *par, const Int_t nn);
  static Double_t Pol1(const Double_t *xx, const Double_t *par);
  static Double_t Pol2(const Double_t *xx, const Double_t *par);
  static Double_t Pol3(const Double_t *xx, const Double_t *par);
  
  static Double_t Landau(const Double_t *xx, const Double_t *par);
  static Double_t InvLandau(const Double_t *xx, const Double_t *par);
  static Double_t ALandau(const Double_t *xx, const Double_t *par);
  static Double_t ALG(const Double_t * xx, const Double_t *par);

  static Double_t RawGaus(const Double_t *xx, const Double_t *par);
  static Double_t ResolutionGaus(const Double_t *xx, const Double_t *par);
  //static Double_t ExpGaus(const Double_t xx[], const Double_t par[]);
  static Double_t TailedGaus0(const Double_t xx, const Double_t amp, const Double_t mean, const Double_t sigma, const Double_t lambda);
  static Double_t TailedGaus(const Double_t xx[], const Double_t par[]);
  static Double_t GausN(const Double_t *xx, const Double_t *par, const Int_t ntype, FFunc basefunc, const Int_t npartype);

  static Double_t HadronBetheBloch(const Double_t *xx, const Double_t *par);
  static Double_t FullBetheBloch(const Double_t *xx, const Double_t *par);
  static Double_t RawAleph(const Double_t *xx, const Double_t *par);
  static Double_t Lund(const Double_t *xx, const Double_t *par);

  /*  
  static Double_t NorExpLandau(Double_t *par);
  static Double_t ExpLandau(const Double_t *xx, Double_t *par);
  static Double_t DimLessExpLandau(const Double_t *xx, const Double_t *pp);
  */

  static Double_t NorDimlessLGE(const Double_t * par);
  static Double_t DimlessALGE(const Double_t * xx, Double_t * par);

  //===================================================================
  //                               Algorithm
  //===================================================================
 public:
  static Double_t GetCorrelatedError(const Double_t v0, const Double_t v1, const Double_t a0, const Double_t a1, const Double_t a2);
  static Double_t GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr);
  static Double_t TruncatedMean(TTree *tt, const TString var, const Double_t lowfrac, const Double_t highfrac, Double_t * grms=0x0, Double_t * gerr=0x0);
  static Double_t TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms=0x0, Double_t * gerr=0x0, Double_t *wws=0x0);

  static void Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint=0);
  static Bool_t FindY(FFunc func, Double_t *par, const Double_t xmin, const Double_t xmax, const Double_t ytarget, Double_t &pos, Int_t &ndepth, const Double_t eps, const Double_t resolution);
  static Bool_t GetMPVFWHM(FFunc func, Double_t *par, const Double_t xmin, const Double_t xmax, Double_t *vout, const Double_t eps=1e-6, Double_t resolution=-999);

  //===================================================================
  //                                Fitting
  //===================================================================
 public:

  enum opt{
    kTMINUIT=0,
    kMINUIT=1,
    kMINUIT2=2
  };

  static TMatrixDSym * GetCov(TMinuit * tmn, const Bool_t ifixed[], const Double_t errs[]);

  static TH1D * GetHfit(const TString hname, FFunc func, const Double_t par[], const Double_t xmin, const Double_t xmax, const Bool_t kbinlogx=0, const Int_t nbins=500);

  static Int_t MaximumLikelihoodFit(const Int_t nx, const Double_t xdata[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t *errs=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);

  static Int_t BinnedLikelihoodFit(const Int_t nx, const Double_t xdata[], const Int_t nbin, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[]=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);
  static Int_t ChisquareFit(const Int_t nx, const Double_t xdata[], const Int_t nbin, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t *errs=0x0, Double_t *chi=0x0, const Double_t thres=0, const Bool_t kxerr=kFALSE, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);

  static Int_t BinnedLikelihoodFit(const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[]=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);
  static Int_t ChisquareFit(const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t *errs=0x0, Double_t *chi=0x0, const Double_t thres=0, const Bool_t kxerr=kFALSE, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);

  static Int_t BinnedLikelihoodFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[]=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);
  static Int_t ChisquareFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t *errs=0x0, Double_t *chi=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);

  static void SetFitPrintLevel(const Int_t lev){ fgFitPrintLevel = lev; }
  static void SetTol(const Double_t tol){fTol = tol;}
  static void SetMaxIter(const Int_t iter){fMaxIter = iter;}
  static void SetMinOpt(const Int_t opt){fMinOpt=opt;}

 private:
  static Int_t FitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);
  static Int_t TMinuitFitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);
  static Int_t MinimizerFitKernel(const std::string minName, const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);

  static Int_t BinnedFit(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Int_t nbin, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);
  static Int_t BinnedFit(const Bool_t kml, const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);
  static Int_t BinnedFit(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix,Double_t *cov, const Double_t *lows, const Double_t *highs);

  static Double_t GetFCNMaximumLikelihood(const Double_t pars[]);
  static void MinFCNMaximumLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);
  static Double_t GetFCNChisquare(const Double_t pars[]);
  static void MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);
  static Double_t GetFCNBinnedLikelihood(const Double_t pars[]);
  static void MinFCNBinnedLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);
  //=================

  static Int_t fgFitPrintLevel;  //print level

  static Int_t fgFitNStep; //n step
  static Int_t fgFitNData; //n data
  static Double_t *fgFitX; //data x
  static Double_t *fgFitY; //data y
  static Double_t *fgFitEx; //err x
  static Double_t *fgFitEy; //err y
  static FFunc fgFitFunc; //fit func
  static Double_t fTol;
  static Int_t fMaxIter;
  static Int_t fMinOpt;
  //===================================================================
  //                               Histograms
  //===================================================================
 public:
  static void PrintHist(const TString tag, const TH1 * hh);
  static void PrintGr(const TString tag, const TGraphAsymmErrors * gr);
  static TGraphAsymmErrors * PublicFilterGraph(const TGraphAsymmErrors * rawgr, const Double_t thres);
  static Bool_t SimpleScaleToRef(const TGraphAsymmErrors * grStat, const TGraphAsymmErrors * grSys, const TGraphAsymmErrors * grRef, TGraphAsymmErrors *& outStat, TGraphAsymmErrors * & outSys, const Bool_t exitOnError=kTRUE);
  static void ScaleGraphPoint(TGraphAsymmErrors * gr, const Int_t ipoint, const Double_t fac);
  static void ScaleGraph(TGraphAsymmErrors * gr, const Double_t fac);
  static void BinLogX(TAxis* ax);
  static TH1D * GetMajorCluster(const TH1D * hin, TH1D * & hrecord, Int_t & clsid, Double_t & cut);
  static Int_t GetNbumps(const TH1D * hh);
  static TGraphAsymmErrors * MeasurementConversion(const TString hn, const TGraphAsymmErrors * gv1, const TH2D * hy2x1, const Bool_t kcutEdge, const Double_t thres);
  static TGraphAsymmErrors * GetPIDCont(TH1D *h1, TH1D *h2);
  static TList * Get2DPIDEff(TH2D *h1, TH2D *h2, const Int_t nxcut, const Double_t xcuts[]);
  static TGraphAsymmErrors * GetPIDEff(TH1D *h1, TH1D *h2);
  //static void GetGrChi2(TList *ll, Double_t grchi2[]);
  static Bool_t GetGrPointLimits(const TGraphAsymmErrors * gr, Double_t & gmin, Double_t & gmax, Bool_t & klogx);
  static Int_t GetCommonLimits(const Int_t ngr,  const TGraphAsymmErrors * grs[], Int_t id0s[], Int_t id1s[]);
  static void GetGrSysMode0(TList *ll, const TString tag, const Bool_t *badgr, const Bool_t kfit, const Double_t hmin, const Double_t hmax);
  static void GetGrSysMode1(TList *ll, const TString tag);
  static TH1D * GetCDF(const TH1D *hh);
  static Double_t TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms=0x0, Double_t * gerr=0x0);

  static void PrintOut(const TString pn, const TCanvas *cc=0, const Int_t kopt=1);
  static void PrintOut(Int_t &id, const TString pn, const TCanvas *cc=0, const Int_t kopt=1);

  static TH2D* NormalHist(const TH2D *h0, const Double_t thres=0, const Bool_t kmax=kFALSE);
  static Double_t * GetAxisArray(TAxis * aa);
  static void FitSlicesY(const TH2D *hh, FFunc ffunc, TH1D *& hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, TH1D *&hchi, const Double_t thres=10, const Double_t lowfrac=0.0, const Double_t highfrac=1.0, TList * ll =0x0);
  static TH1D * GetSeparation(const TH1D *hmpv, const TH1D *hwid, const Int_t anchor, const Bool_t kprint=0, const Double_t gr0=-999);

  static TH1D* ToPDF(const TH1 *h0, const TString hn);

  static Int_t GraphToData(const Int_t ng, TGraphAsymmErrors * grs[], Double_t xmin[], Double_t xmax[],  Double_t xs[], Double_t ys[], Double_t xes[], Double_t yes[]);

  static TGraphAsymmErrors * TreeToGraph(const TString tname, const TString cut, Char_t * sx, Char_t * sy, Char_t * sex0=(Char_t*)"", Char_t * sex1=(Char_t*)"", Char_t * sey0=(Char_t*)"", Char_t * sey1=(Char_t*)"", const Bool_t kXLogToLinear=kFALSE, const Bool_t kprint=kFALSE);

  static void ToGraph(const TString fintxt, const Int_t istart, const Int_t ny, TGraph * gr[]);
  static TGraphAsymmErrors * HistToGraph(const TString grname, const TH1 *hh, const Double_t thres=0, const TH1 *herr=0x0, const Double_t xmin=-1e10, const Double_t xmax=1e10, const Bool_t klogx=kFALSE);

 private:

  //===================================================================
  //                           IO
  //===================================================================
 public:
  static TChain * InputFiles(const TString file, const TString tr, Char_t * dir=0x0);

  //===================================================================
  //                           Physics
  //===================================================================
 public:
  static Double_t ToPt(const Double_t mom, const Double_t eta);
  static Double_t Eta2Y(const Double_t mass, const Double_t pt, const Double_t eta);
  static Double_t DeltaY(const Double_t mass, const Double_t pt, const Double_t eta0, const Double_t eta1);

};


#endif
